﻿void  verletroutine(caja *bx,double Dt)
{
	double posNew[3], dr;
	int index =0;
	double parke=0;
	for(int i=0;i<bx->npart;i++)
	{
		for (int j=0;j<3;j++)
		{
			posNew[j]=(bx->part[i].pos[j])+(bx->part[i].vel[j])*Dt+(bx->part[i].force[j])*Dt*Dt/2.0;
			dr=posNew[j]-(bx->part[i].posOld[j]);
			dr-=int(dr/(0.5*(bx->lbox[j])))*(bx->lbox[j]);
			bx->part[i].vel[j]=(dr)/(2.0*Dt);
			bx->part[i].posOld[j]=bx->part[i].pos[j];
			bx->part[i].pos[j]=posNew[j];
			parke+=(bx->part[i].vel[j])*(bx->part[i].vel[j]);
		}
		bx->boundary(i);
	}
	bx->ekin=0.5*parke;
}

void  verletLinkCellList(caja *bx,double Dt, bool termo)
{
	//double posNew[3];
	double dr;
	double dv[3];
	double parke=0.;
	double posR;
	int nper;
	double anper;
	int mc[3], c;

	///Para calcular la velocidad media del sistema
	dv[0]=dv[1]=dv[2]=0.0;

	//Reset the head of the link-cell list
	for (int i=0;i<bx->ncell3;i++)
	{
		bx->head[i] = -1;
	}

	for(int i=0;i<bx->npart;i++)
	{
		for (int j=0;j<3;j++)
		{
			posR=(bx->part[i].pos[j])+(bx->part[i].vel[j])*Dt+(bx->part[i].force[j])*Dt*Dt/2.0;
			//Condiciones de contorno periodicas
			nper = int(dabs(posR)/(0.5*bx->lbox[j])); //si esta dentro de la caja 0; si esta fuera 1
			anper = sign((double(nper)*bx->lbox[j]),posR); //0 dentro de caja, 1 � -1 si esta fuera
			posR-=anper;

			//posNew[j]=(bx->part[i].pos[j])+(bx->part[i].vel[j])*Dt+(bx->part[i].force[j])*Dt*Dt/2.0;
			//dr=posNew[j]-(bx->part[i].posOld[j]);
			dr=posR-(bx->part[i].posOld[j]);
			dr-=int(dr/(0.5*(bx->lbox[j])))*(bx->lbox[j]);
			bx->part[i].vel[j]=((dr)/(2.0*Dt));
			bx->part[i].posOld[j]=bx->part[i].pos[j];
			//bx->part[i].pos[j]=posNew[j];
			bx->part[i].pos[j]=posR;
		}
		bx->boundary(i);
	}

	///Uso del termostato?
	if(termo)
	{
		vscale2(bx);
	}
	for(int i=0;i<bx->npart;i++)
	{
		for (int j=0;j<3;j++)
		{
			parke+=SQR(bx->part[i].vel[j]);
			dv[j]+=bx->part[i].vel[j];
		}
	

		///Creamos la link-cell list
		for(int j=0;j<3;j++)
		{
			mc[j] = int((bx->part[i].pos[j]+0.5*bx->lbox[0])/bx->lgrid);
		}
		c = int(mc[0]*bx->ncell2+mc[1]*bx->ncell+mc[2]);
		///Link to the previous occupant (or EMPTY if you're the 1st)
		bx->link[i] = bx->head[c];
		///The last one goes to the header"""
		bx->head[c] = i;
		
	}
	for(int i=0;i<3;i++)
		{bx->vel[i]=dv[i]/double(bx->npart);}
	///Introducimos la energia cinetica del sistema
	bx->ekin=0.5*parke/bx->npart;
}

void  LeapFrogLinkCellList(caja *bx,double Dt, bool termo)
{
	//double posNew[3];
	double dr;
	double dv[3];
	double parke=0.;
	double posR;
	int nper;
	double anper;
	int mc[3], c;

	///Para calcular la velocidad media del sistema
	dv[0]=dv[1]=dv[2]=0.0;

	//Reset the head of the link-cell list
	for (int i=0;i<bx->ncell3;i++)
	{
		bx->head[i] = -1;
	}

	for(int i=0;i<bx->npart;i++)
	{
		for (int j=0;j<3;j++)
		{
			posR=(bx->part[i].pos[j])+(bx->part[i].vel[j])*Dt+(bx->part[i].force[j])*Dt*Dt/2.0;
			//Condiciones de contorno periodicas
			nper = int(dabs(posR)/(0.5*bx->lbox[j])); //si esta dentro de la caja 0; si esta fuera 1
			anper = sign((double(nper)*bx->lbox[j]),posR); //0 dentro de caja, 1 � -1 si esta fuera
			posR-=anper;

			//posNew[j]=(bx->part[i].pos[j])+(bx->part[i].vel[j])*Dt+(bx->part[i].force[j])*Dt*Dt/2.0;
			//dr=posNew[j]-(bx->part[i].posOld[j]);
			dr=posR-(bx->part[i].posOld[j]);
			dr-=int(dr/(0.5*(bx->lbox[j])))*(bx->lbox[j]);
			bx->part[i].vel[j]=((dr)/(2.0*Dt));
			bx->part[i].posOld[j]=bx->part[i].pos[j];
			//bx->part[i].pos[j]=posNew[j];
			bx->part[i].pos[j]=posR;
		}
		bx->boundary(i);
	}

	///Uso del termostato?
	if(termo)
	{
		vscale2(bx);
	}
	for(int i=0;i<bx->npart;i++)
	{
		for (int j=0;j<3;j++)
		{
			parke+=SQR(bx->part[i].vel[j]);
			dv[j]+=bx->part[i].vel[j];
		}
	

		///Creamos la link-cell list
		for(int j=0;j<3;j++)
		{
			mc[j] = int((bx->part[i].pos[j]+0.5*bx->lbox[0])/bx->lgrid);
		}
		c = int(mc[0]*bx->ncell2+mc[1]*bx->ncell+mc[2]);
		///Link to the previous occupant (or EMPTY if you're the 1st)
		bx->link[i] = bx->head[c];
		///The last one goes to the header"""
		bx->head[c] = i;
		
	}
	for(int i=0;i<3;i++)
		{bx->vel[i]=dv[i]/double(bx->npart);}
	///Introducimos la energia cinetica del sistema
	bx->ekin=0.5*parke;
}